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Abstract 

Given  an  image  that  has  been  scaled  both  horizontally  and  vertically 
(possibly  with  different  scale  factors  in  the  two  directions),  we 
determine  the  camera  position  and  orientation,  as  well  as  the  scale 
factors  for  sampling  in  the  u  and  v  axes  of  the  image  plane,  from  the 
knowledge  of  the  correspondence  between  a  few  known  points  in  space 
and  their  location  on  the  image  (image  plane).  <  i 
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Abstract 

Given  an  image  that  has  been  scaled  both  horizontally  and  vertically  (  possibly 
with  different  scale  factors  in  the  two  directions ),  determine  the  camera  position  and 
orientation,  as  well  as  the  scale  factors  for  sampling  in  the  u  and  v  axes  of  the  image 
plane,  from  the  knowledge  of  the  correspondences  between  a  few  known  points  in 
space  and  their  locations  on  the  image  ( image  plane ). 

Introduction 

For  many  practical  purposes,  it  is  sufficient  to  model  the  imaging  process  of  a 
camera  by  a  pure  point  projection  process;  then  the  two-dimensional  image  plane 
coordinates  of  a  3-D  point,  as  seen  by  a  camera,  can  be  easily  expressed  in  terms  of  a 
3x4  matrix  using  the  notion  of  the  homogeneous  coordinate  system[2,5].  Because  of 
its  elegance  and  its  usual  accuracy,  this  matrix  has  considerable  applications  in 
fields  like  computer  graphics,  computer  vision,  motion  tracking,  etc.  The  problem  of 
camera  calibration  and  determination  of  the  above  mentioned  matrix  has  received 
considerable  attention  in  the  lilerature[3,4,5,6].  Camera  calibration  is  simply 
determining  the  elements  of  the  matrix  that  describes  the  imaging  process. 
Equivalently,  it  can  be  thought  of  as  determining  certain  geometrical  and  optical 
constants,  such  as  the  direction  of  the  optical  axis,  the  spatial  location  of  the 
camera’s  nodal  point  and  the  focal  length  of  the  lens.  The  solution  that  is  described  in 
what  follows,  is  the  one  suggested  by  Ganapathy  ( [5] ). 


1.  Formulation  of  the  problem. 


Consider  a  right  handed  three-dimensional  orthonormal  coordinate  system 
OXYZ  ( see  Fig.  1 ).  Using  the  homogeneous  coordinate  system  ( [  2  ] )  the  coordinates 
of  a  point  A(  x,  y,  z )  are  represented  by  a  quadruple  of  the  form  (  ax,  ay,  az,  a  )  where 
we  introduced  the  scalar  a  .  In  order  to  obtain  the  three-dimensional  coordinates  of  a 
given  point  that  is  represented  in  the  homogemeous  coordinate  system,  one  merely 
divides  by  the  last  component  of  the  quadruple.  In  a  similar  way,  a  (  two- 
dimensional  )  point  I  (  u,  V  )  of  the  image  plane  (  U-V  )  is  represented  by  the  triple 
(iu,  iv,  i ),  for  some  arbitrary  scalar  quantity  i,  and  we  divide  by  the  third  component 
i  to  obtain  u  and  v. 


Consider  a  camera  having  its  center  at  C  that  is  located  at  {  Xq,  yo,  Zq  )  (  see 
Fig.l);  the  camera  looks  along  a  line  of  sight  (  optical  axis  of  the  camera  )  CO’P 
where  P  is  the  point  at  which  the  optical  axis  of  the  camera  intersects  the  X-Y  plane 
and  O’  the  point  at  which  the  optical  axis  pierces  the  image  plane.  Taking  a  right- 
handed  orthonormal  coordinate  system  0”X’Y’Z’  centered  at  C,  we  can  write  for  any 
point  of  the  3-D  space  the  following  equation: 


—  — 

— 

x’ 

a  b  c 

X 

p 

y’ 

- 

d  e  f 

y 

+ 

q 

z  ’ 

g  h  i 

z 

r 

where 


a  b  c 
d  e  f 
g  h  i 


2 


is  the  rotation  matrix,  and 


P 

q 
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is  the  translation  matrix. 

In  the  above,  (  x,  y,  z  )  and  (  x’,  y’,  z’ )  are  the  coordinates  with  respect  to  the  OXYZ 
and  0”X’Y’Z’  coordinate  systems  respectively. 

Using  homogeneous  coordinates,  the  above  equation  ( 1 )  can  be  written  as  follows: 
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where  we  combined  both  the  translation  and  rotation  matrices  in  one. 

Sincethe(xo,yo,zo)  point  of  the  OXYZ  orthonormal  system  is  mapped  to  the 
( 0, 0, 0 )  point  of  the  O’UVW  system,  then  using  equation  ( 1 )  it  follows  that: 

p  =  -a*xo-b*yo-c*zo  (3) 
q  =  -d*xo- e*yo-f*zo  (4) 
r  -  -g*xo-h*yo-i*zo  (5) 
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U  V:  the  image  plane  coordinate 
system 

OXYZ:  the  world  coordinate  system 
CX’Y’Z’:  the  camera  coordinate 
system 


c*(*o«yo»*o' 


O’ 


Z'.w 


Fig.l 


Substituting  into  equation  ( 1 )  ( x,  y,  z )  by  (  xq,  yo,  zo )  and  ( x’,  y’,  z’ )  by 
(0, 0, 0)  and  solving  for  the  translation  matrix,  we  have  the  following  equation: 


a  b  c 


1  r.i 


But  R,  the  rotation  matrix,  is  assumed  to  be  an  orthonormal  rotation  matrix  and 
therefore  it  follows  that  R"'  =  R'^  .Using  the  latter  equation  we  can  write  (6)  as: 
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This  equation  can  be  equivalently  written  as : 


xo  =  -  p  *  a 

-q*d  -r*g 

(8) 

yo  =  -  P  *  b 

-q*e  -r*h 

(9) 

zo  =  -  p*c 

-q*f  -r*i 

(10) 

Consider  now  a  point  (  x’,  y’,  z’ )  where  the  coordinates  are  taken  with  respect  to 
the  O’X’Y’Z’  system.  The  image  plane  coordinates  of  the  considered  point  will 
therefore  be  given  by  the  equations: 

fl*x’  fl*y’ 

u  =  -  and  V  =  - 

z’  z’ 

Here  we  assumed  that  the  Z’  axis  is  parallel  to  the  optical  axis  and  that  the  image 
plane  is  in  front  of  the  camera;  fl  is  the  focal  length  of  the  camera.  Furthemore, 
X’,Y’,  and  Z’  axes  can  be  aligned  with  the  U,  V,  and  W  axes  such  that  U-axis  is 
parallel  to  the  X’-axis,  V-axis  is  parallel  to  the  Y’-axis,  and  W-axis  is  parallel  to  the 
Z’-axis.  In  practice  though,  the  measurements  in  both  the  X’  and  Y’  axes  are  subject 
to  finite  resolution;  this  consequently  induces  a  discrete  spectrum  of  values  for  the 
measured  quantities.  To  be  general  we  will  assume  two  scale  factors  (possibly 
different)  ku  and  ky  for  sampling  in  the  X’  and  Y’  axes  respectively.  Moreover,  if  (  uq, 


VO,  0  )  are  the  coordinates,  with  respect  to  the  UVW  coordinate  system,  of  the  point 
that  we  consider  to  be  the  one  through  which  the  optical  axis  passes,  then  the  above 
equations  can  be  written  as: 


k^*  fl  *  x’  k^.*  fl  *  y’ 

U  =  Ug  +  — - - ;; -  (11 )  and  v  =  - ; -  ( 12  ) 

Representing  now  the  products  ku  *  fl  and  ky  *  fl  by  ki  and  ka  respectively,  we  can 
clearly  write  the  above  equations,  using  homogeneous  coordinates,  as  follows: 
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(14) 


Note  that  Mi  is  a  3x3  matrix  whose  determinant  is  equal  to  ki  *  ka  ^  0.  Therefore 
it  is  invertible. 
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Observe  that  M2  is  completely  determined  by  the  camera  location  and  orientation. 
Equations  ( 14  )  and  ( 15  )  can  be  written  as : 


where  M3  =  Mi  •M2  is  a  3x4  matrix  that  is  known  as  the  calibration  matrix.  Let 
us  denote  its  components  by  ti  /i  =  1,2,3,4, ...  ,12  .where: 

ti  =  ki*a  +  uo*g 

t2  =  ki*b+  uo*h 

t3  =  ki*  c  +  uo  *  i 


t4  =  ki*p  +  uo*r 


t5  =  k2*  d  +  VO  *  g 


equations  ( 17 ) 


te  =  k2*  e  +  VO  *  h 

t?  =  k2*  f  +  VO  *  i 

ts  =  k2*  q  +  VO  *  r 

t9  =  g 

tlO  =  h 

til  =  i 

tl2  =  r 

Since  we  use  the  homogeneous  coordinates  formulation,  no  change  will  occur  if  we 
divide  all  the  ti’s  by  ti2.  Using  the  notation  Ti  =ti  /  ti2,  we  have: 

Ti  =  (  ki*a  +  uo*  g)/r 
T2  =  (  ki*  b+ uo  *  h  ) /r 
T3  =  (  ki*  c  +  Uo  *  i  )/r 
T4  =  ( ki*  p  +  uo  *  r )  /  r 

T5  =  (  k2*  d  +  Vo  *  g )  /  r  equations  ( 18  ) 

Te  =  (  k2*  e  +  VO  *  h  )  /  r 
T7  =  (  k2*  f  +  VO  *  i )  /  r 


Tg  =  (  k2*  q  +  Vo  *  r  )/r 

Tg  =  g / r 

Tio  =  h  /r 

Til  =  i/r 

Ti2  =  1 

Therefore  equation  (  16  )  caii  be  expressed  in  terms  of  the  eleven  Ti ’s.  Each 
observation  (  u,  v  )  together  with  the  corresponding  (  x,  y,  z  )  gives  us  two  linear 
equations  in  the  eleven  Ti ’s.  Hence  a  minimum  of  six  observations  is  needed  to 
determine  the  entries  of  Mg. 

2.  Some  useful  properties. 

Consider  again  matrix  R  (  rotation  matrix  ).  Since  the  latter  is  an  orthonormal 
rotation  matrix  the  following  identities  hold: 

a^  +  b^  +  c-  =  1  { 19 ) 

d2  +  e2  +  f2  =  1  (  20  ) 

g^  +  h^  +  i^  =  1  (  21 ) 


d*g4-e*h  +  f*i  =  0  (22) 


a*d  +  b*e  +  c*f=0  (23) 


a*g  +  b*h-l-c*i  =  0  (24) 
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Equations  (19)- (21)  represent  the  normality  of  the  vectors  whereas  equations 
(  22  )  -  (  24  )  represent  the  orthogonality  of  the  vectors  that  correspond  to  the  rows  of 
the  rotation  matrix  ( [7] ).  Recalling  that  Det  (  R  )  =  Det  (  R"^ )  =  1  and  that  R'*  = 
R*-,  we  have  the  following  set  of  equations ; 


a  =  e  *  i  -  f  *  h 

(25) 

b  =  f  *  g  -  d  *  i 

(26) 

c  =  d*h-e*g 

(27) 

d  =  c  *  h  -  b  *  i 

(28) 

e  =  c*g  -  a*i 

(29) 

f  =  h*g  -  a*h 

(30) 

g=  b*  f  -  e*c 

(31) 

h=c*d-a*f 

(32) 

i  =  a  *  e  -  b  *  d 

(33) 

3.  Solving  the  problem 

The  equations  (  19  )  through  (  33  )  we  described  above  are  not  all  independent. 
As  a  matter  of  fact  the  existing  interdependence  can  be  expressed  by  the  following 
identity: 

(a*d+  b*e+  c*f)2  =  (  a^  +  b^-i-c ^  )*( d^  +  e^  +  f^  )-(a*e-b*d  )2-(c*d-a*f)2-(b*f-e*c)2 


Consider  now  the  vectors  A,  B,  and  C  where  A,  B,  and  C  are  defined  as  follows: 

A^(Ti,T2.T3)  (34) 


B^(T5.T6,T7) 


(35) 


C  =  (T9,Tio,Tii)  (36) 


Using  the  definitions  of  the  Ti’s,  and  after  some  easy  computations  we  have  that: 


A  A  = 


k2  + 

1  o 


B  B 


,  2  .  2 
•^2  +  ''o 


c  c  = 


A  C  = 


B  C  = 


The  above  equations  clearly  indicate  the  method  that  is  to  be  used  in  order  to 
solve  the  problem. 

3  a)  Computing  the  parameters  of  the  problem 

We  first  determine  the  entries  Ti,  i  =  1,2, 3,4...  11  (  recall  that  Ti 2  =1 ). 

Then  using  the  definitions  (  34  ),  (  35  ),  (  36  )  the  equations  (  37  )  thru  (  41  )  and 
equations  ( 18  ),  we  can  clearly  see  that: 


=  "7~1 — i - T"  ^^^2) 


11 


f  = 


k 


2 


p  = 


(55) 


(56) 


q  = 


(57) 


In  the  above  the  sign  of  r  is  assumed  to  be  positive.  As  a  matter  of  fact  the 
solutions  with  +r  and  -r  correspond  to  the  image  in  front  of  and  behind  the  camera 
respectively  ( see  [5] ).  Moreover,  ki  can  always  be  assumed  positive.  This  is  not  true 
for  kj  though.  The  sign  of  k2  depends  on  the  polarity  of  axes  U  and  V  (  see  [5]  ).  We 
can  determine  the  sign  of  k^  using  the  set  of  equations  { (  28  ),  ( 29  ),  (  30  ) }  and 
{ (  53  ),  (  54  ),  (  55  ) } .  Observe  from  the  latter  set  of  equations  that  the  signs  of  d,  e 
and  f  depend  on  the  sign  of  k2;  therefore,  the  former  set  of  equations  will  hold  only 
for  one  of  the  signs  of  k2 ,  that  is  for  the  correct  one. 

3  b)  How  to  determine  the  table  entries 

Observe  that  we  deal  with  a  total  of  11  unknown  table  entries  and  therefore,  in 
principle,  a  set  of  at  least  51/2  points  is  needed.  Since,  in  general,  more  than  6  points 
are  used,  one  needs  to  apply  the  method  of  least-squares.  Recalling  the  fact  that  we 
are  using  the  homogeneous  coordinate  representation,  we  can  clearly  see  that  the  u 
and  V  coordinates  are  given  by  the  following  equations: 


'  T*x+T  *v+T  •z.+T 
S  I  10  11  1 


u,  = 


12 


T  •  X  +  T^*y  +  T,*z  +  T_ 
5  1  6  ■'i  7  1  8 


V.  = 


’  T*x-l-T  ♦v+T  *z+T 
‘9  *1  ‘lO  *11  ^  ^  ‘12 


/i  =  l,2,3...a  (58) 


/j  =  l,2,3...n  (59) 
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where  s,  w  were  assumed  to  be  equal  to  1  (  see  equation  16  )  and  n  is  the  number  of 
the  experimental  points.  Therefore,  application  of  the  least-squares  method  to  this 
case  will  require  that  we  minimize  the  quantity  A,  where: 


A-X  l(Vx_+Vy,  +  T„-.,f  n'a,-(T,*x,  +T,' y,  tT,*  x +T,  I  + 


I  ((Vx,  +T,„-y,+T„-x+l)'v,-(Vx,+Vy,+T,-x+T,ir 


Taking  the  partial  derivatives  of  A  with  respect  to  Tj,  fori  =  1,2,3,...  11,  and 
equating  them  to  0,  we  come  up  with  the  followingsetof  linear  equations  : 
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T  *(y  x^*u.  )  +  !*(>  x*v.*u  )  +  T  *(y  x.*z.*u.)  +  T/(y  x.*u.)  + 
1  1  1  2^1*11  3  ^  111  4^11 


+  V  ^  S  \  *^i  >  +'^6*  ^  5!  ’‘i*>i*  ''i  V ( X  V ( S  ’‘i*  ^  - 


T  *(y  x2*(u2+v2))-T,  *{  V  X  *y  •(u^+v^))-?,  *(y  X.*z.*(u2+v2))=  Y  x.*(u2+v2)(El) 
9  1  11  10  —  1  ^  1  i  1  11  ^  111  1  ^  11  1 


T  *(  Y  x*y  ♦u  )  +  T  *(  Y  yf*u.)  +  T  *(  Y  y.*z  *  u.  )  +  T/ (  Y  y.*u.)  + 

1  111  Z  11  «}  111  4  1  1 


+  T  *  (  Y  x.*y.*  V  )  +T  *  (  Y  y^*  V.  )  +  T  *  (  Y  y.*z  *  v.)  +  T  *  (  Y  y.*  v. )  - 

5  ^  1  i  1  6  *^1  1  7  ^  1  i  i  8  ^  I  i 


V  ^  X  >>-'^10*'  S  *  (“f  +  vf )  )-Tj,*(  X  y*^*  (uf  +  vf ) )  =  Y  y_*  (u2+  vf )  (E2) 

i=l  1=1  i=l  i=l 


T ,*(  Y  X  ♦z  *u  )  +  T  *(  >'  y  *z.*u  )  +  T  *(  Y  z^*u  )  +  T/  (  z  *  u  )  + 

1  111  2  *^111  3  1  1  4  i  1 


+  T  *  (  Y  X  *z  *  V  )  +T  *  (  Y  y*z  *  V.  )  +  T  *  (  Y  z^*  v.)+T  *  (  Y  z  *  v.  )  - 

5  ^  111  6  ^  *^111  7  ^  j  I  8^11 


T  *(  Y  x.*z.*(u^+v^))-T  *(  Y  y  •z.*(u^+v^))-T,,*(  Y  x^*(u^  +  v^)  =  Y  z  *  (u^+ v^jCES  ) 

9  ^  1111  10  ^  i  1  1  1  11  ^  1  i  1  ^  i  1  1 


T,*(Y  x^)  +  T  ‘(Y  x.*v.)  +  T  *(Y  x*z.)  +  T/(  Y  X.)- 
1-^1  2  ^  i  '  i  3  ^11  4  i 
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-  v>  I  * “i  '  S  “i I  i  “i  I = i  I '  I E-i  I 


V<S",’)'i>  + V'S  yfi  +  V'tyi^'  +  vi  y,'- 


T,‘(l  x,*^i)+V(iy,'z,)  +  T3'(I.f)  +V(i  z)- 


-T,'(S  z,-z'u,I-T,3-(X  y,V“i>-T,2*li  Ilx.-u,!  (E6 


T|''S  X,  l  +  T,-(X  y,>  +  T3MX  z)+T,‘t.- 


n  n  n  n 

-'ra*(  X  X.*u  )-T,-*(  y.*u  )-T,  •(  y  z,*u  )=  y  u 
9  .^11  10  zi-.-'ii  11  ^11 

i=l  i=l  i=l  i=l 


T .♦(y  x^)  +  T .♦(y  X  •y,)  +  T  ♦(Y  x,*z.)  +  T  *(y  X  )- 
5  ^  i  6  I  i  7  i  I  8  ^  i 


-V  ‘  I  1 1  »,  |-T„’<  s  x^-x-  I  =  I  <  X,*  »,  I  I  E9  ) 


■T»*(y  X  ‘y  *v  )  +  T,  *(  y  y^’v  )  +  T,  *{y  y,*z*v  )=  y  (y  *v  )  (E9 

9  ^  i  •' i  \  10  "'ll  11  i  i  i  •' i  I 

i=l  1=1  i=l  i=l 


4.  Estimating  the  '*  quality  ’’  of  the  results 


We  would  like  to  lump  somehow  the  non-linear  nature  ( cf.  equations  ( 17  ) )  of  the 
original  problem  in  a  single  parameter.  At  this  point  one  should  also  recall  the  error 


W'  =  W 


Fig.  2 


introduced  by  the  observations  ( i.e.  (uj,  Vi )  &  (  Xj,  yj,  Zi ) )  made  as  well  as  the  fact 
that  the  camera  used  is  not  a  pin  hole  one  without  distortions.  Assume  that  all  these 
errors  that  are  intrinsic  to  the  experimental  approach  amount  to  the  determination, 
using  equations  (  47  )  -  (  54  ),  to  a  coordinate  system  UVW  such  that  axes  U  and  V 
are  not  perpendicular  to  each  other  but  slightly  off  by  a  small  angle  8  (  see  the  above 
Fig.  2  ).  We  furthermore  considered  axes  U’V’W’,  that  constitute  an  orthonormal 


coordinate  system,  such  that  the  pairs  (  W  ,  W’ )  and  (V,  V’ )  of  axes  represent  axes 
parallel  to  each  other.  It  is  clear  that  the  coordinates  u,v,w  can  still  be  expressed  in 
terms  of  x,y,z  ( see  also  Fig.  1 )  as  follows: 


—  — 

Ml  . 

“  * 

u’ 

a’  b’  c’ 

X 

v’ 

d’  e’  f 

y 

w’ 

g’  h’  i’ 

z 

Ml 

(71) 


where  a’,  b’,  c’,  d’,  e’,  f,  g’,  h’,  i’  are  the  coefficients  of  a  proper  rotation  matrix. 
Moreover  we  observe  that  we  can  write: 


—  — 

M  Mi 

u 

1  0  0 

u’ 

V 

2 

sin5  cos6  0 

v’ 

w 

0  0  1 

w’ 

Mi 

(72) 


Hence  we  can  write  the  following: 


b’ 


b'  •sin6  +  e’  *cos5 


c’ 

c’  *8106  +  1’  ‘cosS 


Comparing  now  equations  (72)  and  (  73 )  we  observe  that  we  can  write: 


a  =  a  ’ 
b  =  b’ 
c  =  c’ 

d  =  a’*sin8  +  d’*cos8 
e  =  b’*sin8  +  e’*cos8 
f  =  c’  *sin8+  f  ’  *cos8 

g  =  g’ 
h  =  h’ 
i  =  i  ’ 

Therefore, 

a*d  +  b*e+c*f  =  sin8*(  a’2  +  b’2  +  c’2 )  +  cos8*(  a’  *  d’  +  b’  *  e’  +  c’  *  f  )  =  sin8 
Recall  that  from  equation  (  23  )  -  "normality”  equation  -  we  have  : 

a*d  +  b*e  +  c*f=0 

Hence,  we  can  estimate  the  "  goodness  ”  of  the  calibration  using  the  value  of 
parameter  8.  It  has  been  observed,  that  the  closer  8’s  value  is  to  0  the  better  the 
calibration  under  consideration.  In  fact,  values  of  8  -  10'^  radians  are  considered 
satisfactory  for  solid  state  as  well  as  vidicon  cameras.  In  case  8  is  significantly 
greater  than  10  it  is  better  to  consider  a  new  set  of  points  than  to  attempt  to 
recalibrate  the  camera  with  the  same  set  of  points. 


5.  Deriving  more  important  parameters  of  the  problem 

So  far  we  have  solved  the  linear  system  that  resulted  after  the  consideration  of 
the  initial  equations  of  the  problem.  In  other  words  we  have  computed  the  entries  of 
the  matrix  M3,  or  equivalently  the  unknowns  Ti’s  /  i  =  1,2,3  ...  12,  (  recall  that 
Ti2  =  1  1  We  still  need  to  translate  these  coefficients  into  camera  parameters, 
namely  the  position  of  the  center  of  the  camera  in  3-D,  the  scale  factors  ku  and  ky  as 
well  as  the  pair  (  uq,  vq  )  ( see  page  6  ). 

The  Uq  and  vq  are  computed  using  equations  (  44  )  and  (  45  ).  Recalling  the 
definitions  ofku  and  ky  we  observe  that  using  the  results  of  equations  (  45  )  and 
(46  )  we  can  compute  ku  and  ky  as  follows: 

ku  =ki  /  n 
ky  =  k2/  fl 

where  fl  is  the  (  known  )  focal  length  of  the  camera. 

Furthermore,  in  what  concerns  the  position  of  the  camera  in  3-D,  one  need  only 
recall  equations  (  8  ),  (  9  ),  and  (  10  ).  These  equations  directly  give  us  the  3-D 
coordinates  of  the  center  of  the  camera. 

At  this  point  it  is  worth  recalling  the  discussion  about  the  sign  of  k2.  As  we  saw, 
according  to  the  set  of  points  considered,  the  sign  of  k2  might  be  different  each  time. 
And  since  the  values  of  d,  e,  f,  and  q  all  depend  on  the  values  of  k2  it  seems  that  we 
must  determine  each  time  the  sign  of  k2,  in  order  to  compute  xq,  yo,  zo-  However  we 
can  observe  that  all  these  variables  appear  in  equations  (  8  ) ,  ( 9  ) ,  ( 10  ) ,  as  products 
of  pairs  and  therefore  their  signs  always  reduce  to  "  -1-  ”.  Hence  there  is  no  need  at 
all  to  compute  the  sign  of  k2. 

Finally,  we  could  also  determine  the  directions  of  the  axes  U,  V,  and  W  with 
respect  to  the  fixed  3-D  coordinate  sytem  OXYZ.  In  order  to  do  this  we  only  need  to 
consider  the  (  orthonormal )  vectors: 

ni  =(a,b,c) 


na  =(g.  h.i) 

As  a  matter  of  fact  ni  is  the  vector  that  is  parallel  to  U-axis,  n2  is  the  vector  that 
is  parallel  to  the  V-axis  and  na  is  the  vector  that  is  parallel  to  the  W-axis. 

Note:  As  it  has  been  also  pointed  out  by  Ganapathy  ( [5]  ),  the  computation  ofuo.vo  is 
rather  sensitive.  Therefore,  if  one  is  interested  in  a  quite  precise  computation  of  those 
quantities,  one  should  follow  another  approach  ( see  APPENDIX  I ). 

6.  Deficiencies  of  the  algorithm 


The  above  described  algorithm  has  one  important  deficiency.  Assume  that  the 
considered  set  of  points  has  a  symmetry  with  respect  to  any  of  the  planes  defined  by 
any  two  axes  of  the  3-D  coordinate  system.  If  the  optical  axis  lies  on  this  plane,  then 
some  of  the  summations  appearing  in  equations  (  El )  through  ( Ell )  will  be  equal  to 
0.  Which  ones  exactly  are  equal  to  0,  depends  on  the  symmetry.  However,  this  is 
enough  to  reduce  the  rank  of  the  coefficients’  matrix  (  matrix  Aij  below  ).  But 
reduction  of  the  rank  of  this  matrix  has  as  a  consequence  our  inability  to  solve  the 
linear  system  of  the  11  equations  in  the  11  unknowns  Tj,  i  =  l,  2,  3,  ,.  11.  It  is 
therefore  clear  that  one  should  not  choose  such  configurations  of  points.  However, 
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the  problem  is  not  necessarily  solved.  This  fact  lies  on  more  intrinsic  reasons.  Let  us 
consider  the  above  linear  system  of  equations. 
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Definition:  The  above  system  of  equations  is  called  ill-conditioned  if  the 
determinant  of  the  matrix  of  coefTicients  has  a  value  that  is  very  close  to  0. 


In  such  cases,  an  error  is  introduced  in  the  values  of  the  computed  parameters. 
Moreover,  consider  the  case  where  the  coefTicients  Ay  are  measured  with  an  error 
±AAij.  This  means  that  the  value  of  Ay  lies  in  an  interval  [A*y-ey  ,  A*ij  +ey], 
/  i  =  1,2,3, ...  n  &  j  =  1,2,3, ...  m.  The  "stared”  values  A*ij  are  the  ones  that  "correct” 
measurements  should  give.  It  is  clear  that  we  cannot  be  sure  that  for  no  set  of  values 
Ay,  the  above  system  become  ill-conditioned. 

Definition:  We  call  a  linear  system  critically  ill-conditioned  if  it  is  not  ill- 
conditioned,  but  there  is  at  least  one  set  of  values  for  Ay  in  the  intervals  [  A*y  -Cy, 
A*ij  -f-Cy  ],  /  i  =  1,2,3,  ...  n  &  j  =  1,2,3,  ...  m,  for  which  the  linear  system  becomes  ill- 
conditioned. 

In  a  recent  work  ( [1] ),  that  dealt  with  this  problem,  a  measure  of  ill-conditioning 
was  determined.  In  fact,  the  quantity: 

n  m 

(b  =  ^  ^  I  b  .  I  •  E  .. 

i  =  1  j=  1 

where  by  =  (  (  A'l  )^  >  ij  /  i  =  1,2,3,  ...  n  &  j  =  1,2,3,  ...  m,  can  be  used  to  determine 
whether  the  considered  system  is  ill-conditioned.  The  following  hold: 

(J)  <  1  =>  Not  critically  ill-conditioned 
<J)  >  1  =>  Critically  ill-conditioned 

Furthermore,  an  upper  bound  for  the  error  AXi  /  i  =  1,2,3  ...  m,  in  the  computed 
values  of  Xk/k  =  1,2,3, ...  m,  was  established; 


where  ||  a  ||  i  is  the  first  norm  of  matrix  a  =  (  ajk  )  /  j  =  1,2,3, ...  n  &  k  =  1,2,3, ...  m,  that 
is  defined  as  follows: 

n  m 

II  a  II  j  ^  I  a.^  I  }or  max^  {  ^  I  a I } 

j=l  k=l 

One  could  use  the  parameter  to  determine  whether  the  system  of  equations 
( El )  thru  (Ell  ),  is  ill-conditioned  or  not,  and  therefore  determine  the  "  quality  ”  of 
the  calibration.  However,  the  necessary  overhead  does  not  make  the  effort  worth  it, 
since  one  can  easily  tell  from  the  calibration  matrix  whether  the  system  is  ill- 
conditioned  ( in  which  case  most  of  its  entries  will  be  equal  to  0  ),  or  not. 

7.  Implementation  of  the  described  algorithm 

The  algorithm  we  just  descibed,  has  been  implemented  in  the  University  of 
Rochester,  during  the  spring  of  1986.  A  full  description  of  the  program’s  output 
together  with  useful  hints  and  other  relevant  information  can  be  found  in 
APPENDIX  II. 
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APPENDIX  I 


The  following  describes  a  quick  method  for  the  determination  of  the  of  the 
coordinates  of  the  pixel  on  the  DataCube  display  through  which  the  optical  axis  of 
the  camera  in  use  passes.  The  idea  of  the  method  relies  on  the  notion  and  the 
properties  of  the  so-called  'Vanishing  point”. 

Method:  A  set  of  8  points  forming  the  vertices  of  a  cube  is  constructed  in  3-D  space 
(see  Fig.  l.a  ).  These  points  can  be  constucted  using  8  nuts  screwed  on  vertical 
threaded  rods,  or  using  a  set  of  8  targets  like  the  one  shown  in  Fig.  1  .b.  The  camera 
under  consideration  is  moved  until  the  points  A,  A’,  B,  B’,  C,  C’,  D,  D’  of  the  image 
plane  (  Datacube  display  )  that  correspond  to  the  "vertices”  of  the  original  cube,  are 
such,  that  the  quadruples  (  A,  A’,  C,  C’ )  and  (  B,  B’  D,  D’ )  represent  sets  of  points 
lying  on  two  distinct  lines  of  the  image,  say  ei  and  e2  respectively  (  see  Fig.  2. a  and 
Fig.  2.b  ).  Lines  ei  and  e2  intersect  at  a  point  P  (  the  vanishing  point  ).  The 
coordinates  of  this  point  are  the  coordinates  we  are  after.  How  can  one  determine  the 
coordinates  of  P  ?  Recall  that  A,  A’,  C,  C’  lie  on  the  same  line,  namely  ci.  Using  a 
cursor  moving  routine,  one  determines  the  coordinates  of  points  A,  A’,  C,  C’.  Then 
using  a  simple  calculator,  and  applying  a  least-squares  fitting  method,  the  slope  and 
the  offset  of  line  ei  are  computed.  The  same  procedure  is  repeated  for  the  case  of  line 
£2.  The  slopes  and  offsets  of  lines  ei  and  £2  being  determined,  the  coordinates  of  P 
can  be  computed;  in  fact  one  need  only  express  £i  and  £2  as  linear  equations  in  two 
unknowns,  say  x  and  y,  and  then  solve  the  resulting  linear  system.  The  solution  of 
this  system  will  be  the  pair  of  coordinates  of  P  . 

Improvement  of  the  results  can  be  achieved  by  repeating  the  above  procedure  using 
other  rigid  configurations  ( see  Fig.  3.a  and  Fig.  3.b ). 

Note.-  Obviously,  the  above  described  method  is  general,  and  can  be  applied  to  any 
camera  and  display  combination. 

Results:  Application  of  the  above  described  method  for  the  case  of  a  cube  gave  the 
following  results,  for  the  cameras  of  our  laboratory; 

CAMERA  TYPE  COORDINATES 

vidicon  (  353,  247  ) 

ccd  ( 334,  200  ) 


APPENDIX  II 


This  appendix  gives  a)  some  very  useful  observations  that  have  been  made  as 
well  as  hints,  and  b)  the  file  structures  that  are  needed  together  with  a  description  of 
the  output  of  the  actual  program. 

Part  I 

Choosing  the  3-D  points 

This  is  the  most  crucial  part  of  the  whole  procedure.  In  fact  the  selection  of  the 
world  points  can  lead  to  a  disaster  or  a  success.  Some  useful  observations  have  been 
made: 

a)  choose  a  sufficient  number  of  points;  although  the  algorithm  requires  only  5  1/2 
points,  a  "trustworthy”  calibration  was  found  to  require  at  least  8  or  9  points. 

b)  pay  attention  to  the  configuration  of  the  3-D  points  with  respect  to  the  optical  axis 
and  the  3-D  (  fixed  )  coordinate  system;  besides  some  obvious  degenerate  cases  for 
which  the  algorithm  will  not  work,  there  are  also  other  cases  that  may  cause 
problems.  As  a  general  suggestion,  though,  one  should  keep  in  mind  that  the  more 
randomly  the  targets  are  selected,  the  higher  the  possibility  for  a  succesful 
calibration.  It  is  obvious,  of  course  that  one  should  not  select  configurations  that 
experience  some  kind  of  symmetry  (  especially  axial-symmetry  )  with  respect  to  the 
optical  axis. 

c)  the  3-D  points  should  be  selected  "a.s  far  as  possible”  from  the  camera  under 
consideration;  of  course  this  is  in  direct  relation  with  the  focal  length  of  the  camera. 
In  general,  the  targets  must  be  placed  at  such  positions,  so  that  all  of  them  are 
focused.  Obviously,  if  the  camera  supports  apperture  setting  facilities,  one  would 
like  to  set  the  apperture  to  16  or  22  under  intense  lighting  conditions. 

Measuring  coordinates 

®  3-D  space:  It  has  been  observed  that  measuring  the  3-D  coordinates  of  the 
world  points  with  an  accuracy  of  approximately  0.1  inches  is  sufficiently 
satisfactory.  However,  one  should  keep  in  mind  that  all  the  quantities  that  have  a 
dimension  of  MOL^TO,  e.g.  focal  length  of  the  camera,  world-point  coordinates  etc., 
must  be  expressed  in  the  same  length  units,  e.g.  all  in  inches  or  all  in  millimeters ... 

®  Image  plane:  Using  a  cursor  moving  routine,  one  must  determine  the  image 
plane  coordinates  with  an  inaccuracy  of  at  most  ±3  pixels. 


Discovering  potential  deficiencies  of  the  camera  in  use 

One  problem  that  might  arise  is  the  following:  "cheap  optics”  or  "not  properly 
mounted  on  the  camera”  lens.  One  could  easily  figure  out  whether  this  is  a  problem 
for  the  camera  in  use,  or  not,  by  placing  some  targets  in  front  of  the  camera  and 
trying  to  focus  on  them;  if  their  images  are  displaced  in  a  non-consistent  way,  then 
there  is  a  high  possibility  of  the  camera  having  the  above  mentioned  problem.  Of 
course  it  is  needless  to  say  that  there  is  no  hope  of  calibrating  this  camera  . 

Algorithm’s  results 

It  has  been  observed  that  the  computation  of  the  values  uq,  vq,  is  rather  sensible 
and  related  to  the  direction  of  the  optical  axis.  In  fact,  one  should  use  another 
approach  in  order  to  get  more  robust  values  for  these  parameters  ( see  APPENDIX  I). 
Besides  that,  the  remaining  parameters  are  precisely  computed.  It  must  be  pointed 
out,  though,  that  if  one  wants  to  get  reliable  results,  one  must  try  a  few  calibrations 
with  the  camera  being  horizontal,  and  for  different  camera  positions  (  with  respect  to 
the  3-D  coordinate  frame  )  and  sets  of  world-points,  in  order  to  get  a  feeling  about  the 
the  response  of  the  camera  in  use.  Then  other  positions  might  be  tried  (  pan,  tilt  or 
roll  3;  0 )  as  well. 

Part  II 

Description  of  the  necessary  file  structures 

The  files  that  are  used  by  the  calibration  routine  are  two: 
i)  the  file  "  data3D  ”,  that  contains  the  number  of  experimental  points  as  well  as 
their  3D  coordiates.  This  file  must  be  set  up  by  the  user,  and  should  have  the 
following  format: 

line  #1  #of _ points _ used 

line  #2  Xi  ( space )  Yi  ( space  )  Zi 
line  #3  X2  ( space )  Y2  ( space  )  Z2 
line  #4  X3  ( space )  Y3  ( space )  Z3 

etc. 

Note:  The  Xi’s,  Y  i’s  and  Zi's  are  read  as"  double’s 
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ii)  the  file  "  dataimage  ”  that  contains  the  image  plane  coordinates  of  tht 
experimental  points.  This  file  is  either  set  up  by  the  user,  in  which  case  it  has  u» 
have  the  following  format; 

line  #1  ui  (  space  )  vj 
line  #2  U2(  space  )V2 

line  #3  U3  (  space  )  V3 
line  #4  U4(  space  )v4 


or  by  invoking  the  special  command  called  "  dqstorepoint  the  latter  is  a  cursor 
moving  routine  that  uses  the  "mouse”  of  a  SUN-Workstation  that  is  connected  to  the 
DataCube.  The  user: 

1)  selects  points  by  clicking  the  left  button  of  the  mouse 

2)  confirms  a  selected  point  and  makes  an  entry  for  the  corresponding  coordinates  to 
the  file  "  dataimage  ”,  by  clicking  the  middle  button  of  the  "mouse" 

3)  "  exits  ”  by  clicking  the  right  button  of  the  "mouse”. 

Note:  It  is  worth  mentioning  that  the  entries  of  the  ^  i-h  1  )-th  line  of  the  file  "dataM) 
and  the  ones  of  the  i-th  line  of  the  file  "dataimage"  correspond  to  the  i  th  expennienta! 
point,  i  =  l  2,3,4,5,6, ...  n. 

A  typical  session 

The  following  is  a  script  of  a  typical  session.  We  assume  that  the  files  "  data31)  ' 
and  "  dataimage  ”  were  appropriately  set  up  as  described  abr  The  names  of  the 
parameters  are  consistent  with  the  ones  appearing  to  the  description  of  the 
algorithm... 

<  ISIDORE  >>>  %  cat  data30 
9 

552.2  376.5  1989.1 
815.75  376.5  1989.1 
815.75  169.4  1989.1 
351.4  169.4  1689.3 
351.4  376.5  1689.3 
702.8  169.4  1185.975 


702.8  376.5  1185.975 
753.0  169.4  1085.575 
753.0  376.5  1085 .575 
<  ISIDORL  >>>  %  cat  dataimage 
286  164 

504  164 

505  318 
6?  32  7 
64  143 
43?  37? 

429  89 
517  385 
5  17  6  7 

'  ISIDORE  >>>  7.  dqcalibrate 


Please  enter  the  focal  length  of  the  camera  in  use:  25 


Reading  the  3-D  coordinates  Number  of  points=  9 


Ru 1 n  t  numte  ^  1 

55?  200000  376.500000  1989.100000 

Point  number  2 

815.750000  376.500000  1989.100000 

Point  number  3 

815.750000  169.400000  1989.100000 

Point  number  4 

351.400000  169.400000  1689.300000 


351.400000  376.500000  1689.300000 

Point  number  6 

702.800000  169.400000  1185.975000 

Point  number  7 

702.800000  376.500000  1185.975000 

Point  number  8 

753.000000  169.400000  1085.575000 

Point  number  9 

753.000000  376.500000  1085.575000 

Reading  the  image-plane  coordinates 

Point  number  1 
286.000000  164.000000 

Point  number  2 
504.000000  164.000000 

Point  number  3 
505.000000  318.000000 

Point  number  4 
62.000000  327.000000 

Point  number  5 
64.000000  143.000000 

Point  number  6 
432.000000  372.000000 


429.000000  89.000000 


i 

iM 


I*"' 


Point  number  8 
517.000000  385.000000 

Point  number  9 
517.000000  67.000000 


The  solution  matrix  is  the  following  : 

-5.288670  -0.060201  -1.428586  3957.598778 

-0.051500  4.650256  -0.929366  -914.765229 

-0.000063  -0.000189  -0.003655  1.000000 

r2  =  74628.615332 


Observe  ad  +  be  +  cf  =  0.011060 

The  skew  angle  estimation  parameter  is  equal  (in  radians)  to 
;  0.011060 


u0=415. 503321 


v0=188. 044822 


The  resolution  parameters  are: 
Ku  =  57. 512410  Kv  =  51. 273858 


The  remaining  parameters  are  : 

a  =  -0. 999847  b  =  0. 003509  c  =  0. 017124 


d  = -0.008440  e  =  0. 998633  f  = -0.051583 
g  =  -0. 017282  h  =  -0. 051723  i  =  -0.998512 


I 

I 


i 


A 

i 

■i 

I 

A 


: 


p=672. 994285  q=-235 . 026810  r=273. 182385 


The  directions  of  the  3  axes  with  respect  to  the  fixed  3D 
coordinate  system  are 

U-axis:  (-0.999847.0.003509,0.017124) 

V-axis:  (-0.008440,0.998633,-0.051583) 

(optical )Z -axis:  (-0.017282.-0.051723.-0.998512) 


The  position  of  the  center  of  the  camera  is  at: 
Xo=675. 628994  Yo=246 . 473854  Zo=249 . 128475 


<  ISIDORE  »>  % 


